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We analyze the dynamical and energetic instabilities of spin currents in a system of two-component 
bosons in an optical lattice, with a particular focus on the Neel state. We consider both the weakly 
interacting superfluid and the strongly interacting Mott insulating limits as well as the regime near 
the superfluid-insulator transition and establish the criteria for the onset of these instabilities. We 
use Bogoliubov theory to treat the weakly interacting superfluid regime. Near the Mott transition, 
we calculate the stability phase diagram within a variational Gutzwiller wavefunction approach. In 



the deep Mott limit we discuss the emergence of the Heisenberg model and calculate the stability 
diagram within this model. Though the Bogoliubov theory and the Heisenberg model (appropriate 
£N| for deep superfluid and deep Mott phase respectively) predict no dynamical instabilities, we find, 

interestingly, between these two limiting cases there is a regime of dynamical instability. This 
O result is relevant for the ongoing experimental efforts to realize a stable Neel-ordered state in multi- 

component ultracold bosons. 
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I. INTRODUCTION 

Ultracold atomics gases have recently emerged as a 
very important platform to study non-equilibrium quan- 
tum dynamics of interacting many-body systems. The 
tunability of Hamiltonian parameters together with al- 
most complete isolation from the environment and the 
long time-scales in these systems have made it possible 
to study intrinsic non-equilibrium dynamics of these sys- 
tems without ultrafast probes. Recently, there has been 
a growing number of experimental and theoretical inves- 
tigations of the dynamical properties of Bose-Einstein 
Condensates (BECs) in optical lattices [1]. Of partic- 
ular interest are experiments which exhibit a dynamical 
instability, which is a generic phenomenon present in non- 
linear systems under appropriate conditions. Previously 
established examples of dynamical instabilities occur in 
water waves [2, 3], light in dielectric media [4-7], and 
plasmas [8-10]. Recently, dynamical instabilities have 
been observed in ultracold fermi gases after a tuning of 
the interaction parameters [11, 12]. The realization of dy- 
namical instabilities for current-carrying states in BECs 
have received considerable attention both theoretically 
[13-19] as well as experimentally [20-27]. 

Two qualitatively distinct types of instabilities can oc- 
cur for interacting systems of bosons: (i) the energetic 
instability and (ii) the dynamical instability. The en- 
ergetic instability occurs if the system is not at a local 
minimum of the mean-field energy. If the system is ca- 
pable of dissipating energy, then it will decay from the 
initial metastable state thereby exhibiting the instabil- 
ity. A well-known example for this case is the Landau 
instability (LI) for which a superfluid carrying current in 
excess of the sound velocity becomes unstable, leading 
to a breakdown of superfluidity. In contrast, a dynami- 
cal instability (DI) occurs when the system has collective 
modes with complex frequencies. Such modes will result 
in an exponential growth of small perturbations which 



manifests as a rapid depletion and fragmentation of the 
condensate [22, 27]. For systems that do not have a dis- 
sipative mechanism, the energetic instability alone will 
not occur. On the other hand, the dynamical instabil- 
ity occurs even without dissipation, and will be observ- 
able unless the growth time of the most unstable mode is 
longer than experimental time scales. It can also be seen 
that an energetic instability is a necessary condition for 
a dynamical instability. 

Bosons in an optical lattice undergo a quantum phase 
transition from a superfluid phase in the weakly inter- 
acting limit to an incompressible Mott insulator phase 
as the interaction parameter increases beyond a critical 
value. Scalar bosonic condensates have a U(l) symme- 
try associated with the superfluid phase, resulting in a 
conserved mass current. In the presence of a lattice, 
when the externally imposed current exceeds a critical 
value, the system manifests a dynamical instability [13- 
15]. The critical current required for the dynamical in- 
stability decreases with increasing interaction strength 
and vanishes at the critical interaction required for the 
superfluid-insulator transition. Additional types of dy- 
namical instabilities can occur in multicomponent con- 
densates due the their more complex order parameters. 
In this paper we will focus on two-component bosons in 
an optical lattice with spin-independent interactions. In 
addition to the superfluid-insulator transition, this sys- 
tem also shows a spontaneous ferromagnetic spin order- 
ing in the equilibrium ground state. This system has an 
SU(2) symmetry due to invariance of the energy under 
spin rotation which results in a conserved spin current. 
In the presence of externally imposed spin currents (spin 
twists), the system exhibits dynamical instabilities when 
the spin current exceeds a critical value. We will mainly 
focus on these spin-current driven instabilities, which oc- 
cur in addition to and even in the absence of any mass 
current driven instabilities. 

Previous work addressing spin current instabilities in 
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bosonic systems have focused on the continuum, weakly 
interacting superfluids where the Gross-Pitaevskii Equa- 
tion (GPE) is applicable. In such a context, the coun- 
terflow instability [28-32] as well as the instability of a 
spin-one condensate from an initial helical state [33, 34] 
have been investigated. Here we analyze the instabili- 
ties of the system in the presence of an optical lattice 
for a wide range of interaction parameters going from 
the weakly interacting limit (the deep superfluid phase) 
through the intermediate regime near the superfluid- 
insulator transition to the strongly interacting (atomic) 
limit, deep into the Mott phase. The weakly interacting 
regime is treated within the standard Bogoliubov theory, 
while the strongly interacting regime is treated within 
a spin-wave approximation of the ferromagnetic Heisen- 
berg model, where the spin-spin interaction comes from 
super-exchange mechanism. The intermediate interac- 
tion regime is treated within a variational Gutzwiller 
wavefunction ansatz. We extend the Gutzwiller ansatz 
to both the deep superfluid and the deep Mott limit and 
compare the results with those from the more established 
formalisms mentioned above. 

To analyze stability of the bosonic states, we con- 
struct either mass or spin current carrying mean-field 
states. The spectrum of quantum fluctuations about 
these stationary states is then calculated within a Gaus- 
sian approximation. Negative eigenvalues of the fluctu- 
ation Hamiltonian indicate an energetic instability while 
a complex collective mode spectrum indicates a dynam- 
ical instability. For a dynamical instability, the positive 
imaginary part of the complex spectrum gives the growth 
rate of the unstable fluctuation modes. Our main results 
are: (i) We show that the mass current induced instabili- 
ties give rise to the same instability phase diagram in the 
critical current interaction plane for both spinless and 
two-component bosons, (ii) The two-component bosons 
exhibit a spin-current induced dynamical instability in a 
large region of the critical current interaction strength 
plane in the superfluid phase. We also show the collec- 
tive modes which are unstable and compute their growth 
rates, (hi) We focus on the Neel-ordered state, which 
can be interpreted as a spin-current-carrying state with 
particular commensurate wave vector. Although the Neel 
configuration is not the ground state of the system, there 
are proposals [35] to experimentally explore the physics 
about this high energy state provided its lifetime is suf- 
ficiently long. We show that while this state is stable 
in the deep superfluid and insulating limits, in the inter- 
mediate regime, interestingly, the system is dynamically 
unstable. We thus provide a comprehensive picture of 
the spin-current induced dynamical instabilities in two- 
component bosons on optical lattices for a wide range of 
interactions and spin currents. 

The paper is organized as follows. In Sec. II, we re- 
view the established DI of the mass current of spinless 
bosons. We use the Bogoliubov theory to analyze the 
superfluid limit and the Gutzwiller ansatz to analyze the 
strongly-interacting regime close to Mott boundary. This 



prepares us to investigate the instabilities related to the 
spin current of a two-component Bosonic condensate in 
Sec. Ill in the regime of weak as well as intermediate in- 
teractions. We shall present the stability phase diagram 
and discuss how our results are connected to the deep 
Mott limit. In Sec. IV we discuss the stability of the 
Neel state limit for different regimes. Finally in Sec. V 
we summarize our results. 



II. INSTABILITIES OF MOVING SCALAR 
CONDENSATES 

For completeness and to set the notation and general 
approach, we first briefly consider the mass current in 
a single-component BEC and the concomitant Landau 
and dynamical instabilities. The weakly interacting su- 
perfluid case was originally considered in Refs. [13, 15], 
while the regime near the Mott transition was addressed 
in Ref. [18, 19]. A system of bosons on a lattice and in 
the lowest band is described by the Bose-Hubbard model 

H = -* E + hc ) +| E ^ - ^-i 1 E ( n » - ») 

(ij) i i 

(i) 

where b\ is the boson creation operator on the lattice 
site i, t is the hopping matrix element, U is the on-site 
repulsion, \i is the chemical potential, ni = b\bi, and n 
is the average number of bosons per site. We consider 
this model in one, two, and three dimensions for cubic 
lattices. When t ^> U/n the system has a superfluid 
ground state, and Bogoliubov theory describes its ele- 
mentary excitations. When t ~ U/n, there is a quantum 
phase transition at U = U c to an incompressible Mott 
state. The Bogoliubov theory fails in the vicinity of this 
transition, however, a variational Gutzwiller ansatz can 
be used to treat the system in this regime. 



A. Weakly Interacting Superfluid 

Deep in the superfluid phase the current-carrying 
states can be represented by a condensate wavefunction 
of the form 

bi = vW p ' x *. (2) 

which has a phase twist along p and carries a mass cur- 
rent oc sin p • Xij between neighboring sites. This wave- 
function can be found within mean-field theory by solving 
the Gross-Pitaevskii equation. Expanding the energy of 
the system, Eq. (1), about this state to quadratic order, 
with bi = bi~\-(j)i, one obtains the fluctuation Hamiltonian 

$H = E k *k M k(p)*k where = (4,4>_ k ) and 
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Phase twist 


Spin Twist 




Continuum 


Lattice 


Continuum 


Lattice 


LI 


p > \Jvn\Jn 


sin 2 p ^ Un 
cosp zt 




p ^ mod 27r 


DI 


Never 


cosp < 


p#0 


pzfz 0,=bf mod 2tt 



Table I: The conditions for Landau and dynamical insta- 
bilities of mass-current carrying states (for spinless bosons) 
and spin-current carrying states (for 2-component spin-full 
bosons) in the weakly interacting limit on a square lattice 
(calculated within Bogoliubov theory). The results are valid 
for one, two, and three dimensions. The case of twisting along 
the diagonal of the square lattice: p = p^2 i Xi is taken. For 
comparison, the conditions for the instabilities in the contin- 
uum are also given. 



with e q = — ztjq, where z is the coordination number and 
7 q = z~ x e* q (5 . The energies of the normal modes of 
the system are given by the eigenvalues of the matrix 
<j z Mk [13] where <j z is a Pauli matrix. On the other 
hand, if the system is at a local minimum in energy, then 
the matrix itself will be positive definite. We thus 
summarize the following criteria for the instabilities: 

• LI: at least one eigenvalue of is negative 

• DI: at least one eigenvalue of a z M\^ is complex 

For mass-current-carrying states, it is well known that 
the continuum theory only sustains Landau instabilities, 
which occur when the current in the system exceeds the 
speed of sound. There are no dynamical instabilities in 
the continuum theory. However, on a lattice the sys- 
tem exhibits both Landau and dynamical instabilities 
with the criterion for critical current summarized in Ta- 
ble I. The dynamical instability is crucially related to the 
softening of collective modes at finite wavevectors, which 
does not occur in the continuum. 



B. Gutzwiller Ansatz 

To investigate the DI for the Bose-Hubbard model 
for stronger interactions, we shall approach the 
problem within a truncated Hilbert space. We con- 
sider the variational Gutzwiller wavefunction for the 
ground state, \Q) = Yl i with |^) = cos ^ \n) i + 
n-lj.+sinf e i V i + 



sin - 



[cos 



where \n) i are the Fock states on the site i. This 
variational state was used in Ref. [36] to study the 
Bose-Hubbard model near the Mott transition in the 
absence of a current. Our calculations follow along 
similar lines with an important distinction: the phase 
ifi is position-dependent, i.e. cpi = p • such that 
(0| bi |0) oc sin#e 2pXi , which ensures that the state 
carries a mass current flowing along p. Other parameters 
are then varied to minimize the energy of this mean-field 
state, giving X = f an d rj = 0. 




0.0 0.2 0.4 0.6 0.8 1.0 

u/u c 

Figure 1: Stability phase diagram of phase twists in a lattice 
in different dimensions, where twisting is imposed along: p = 
px . The same diagram is obtained for either spinless bosons 
Eq. (1) or the two-component bosons Eq. (5). This diagram 
agrees well with the more extensive treatment presented in 
[18, 19]. 



The Hamiltonian is expanded about this stationary 
state in the following way: we introduce the bosonic pseu- 
dospin operators t^ ai |vac) = \n-\- a) v a = { — 1,0, 1} with 
the constraint ^2 a t^ ai t a i = 1, so that the Boson opera- 
tors can be written as b\ — > y/Et^t-u + y/n-\- listen- A 
unitary transformation is then performed with 




-*p-*i . 

o-*p-x» 6 



cos - 



sin i 



COS 7 



\ 



V 



V2. 







/it 



7 



li 



v 0i 

(4) 

and the Hamiltonian is written in terms of the d op- 
erators. Since |vac) represents the minimum energy 
state, it is macroscopically occupied, while d) n>Q are fluc- 
tuations about this state. Therefore, we eliminate d\ 
using d\ « do ~ 1 — ^d\ { dn — \S 2i d2i^ which resembles 
the Holstein-Primakoff transformation [37] used in spin 
models. 

To the quadratic order in the operators 4l>o> the 
Hamiltonian has the form H = ^M^^k where ^ = 
(d\ k1 c?i,-fe? <4fc> <^2,-fc^ and the form of Mk is given in 
Appendix A. For a given U/t and p, we compute the 
energies oJi^k for k G [— 7r,7r] by a Bogoliubov transfor- 
mation. As noted before, the presence of complex eigen- 
frequencies indicate a dynamical instability. 

For direct comparison with previous work, we con- 
sider p along an axis of a d-dimensional cubic lattice 
(p =pxi), giving 7 P = cos;p ^ d ~ 1 ^ ) . The resulting phase 
diagram is shown in Fig. 1, which shows good agreement 
with the results in [18, 19] where a numerical analysis is 
performed, taking a larger Hilbert space. 

The Bogoliubov analysis is justified only if the fluctu- 
ation occupation (d^du + dj^cfei^ is small compared to 
unity. This is checked in the stable regimes after the 
Bogoliubov transformation is done. We find that for 
the 2D system, the fluctuation is less than 0.2 for all 
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U > 0.1U C and reach up to 0.5 as U —> and p —> tt/2. 
This means the quantitative result should be trustworthy 
for U > 0.1U C . However in ID we always find divergent 
occupation of the fluctuations as expected because of the 
significance of quantum fluctuations. The qualitatively 
good agreement for ID results with experiment might 
be understood as due to the logarithmic nature of the 
divergence, which is not severe in finite-sized systems. 



III. SPIN CURRENT INSTABILITIES IN 
TWO-COMPONENT CONDENSATES 



Having set up the formalism to study Landau and dy- 
namical instabilities in spinless bosonic systems, we will 
now adapt this formalism to study instabilities of spin- 
current carrying states in condensates of two-component 
bosons. The starting point for our analysis is the two- 
component rotationally invariant Bose-Hubbard model 



form 



H 



-/i^(n; -n) 



bjcj + he 



U 



( n i - n Y 



(5) 



where b\ a creates a boson of spin a on site z, = 

XI a fy[cAcr> and n is the average particle number per site. 
Such a system could be realized using, for instance, two 
hyperfine states of alkali atoms. Due to the smallness of 
spin-exchange interaction for typical alkali atoms, such 
systems possess an approximate SU(2) symmetry, which 
is reflected in the spin-independent form of the inter- 
actions which we consider here. For simplicity we will 
concentrate on the case when n = 1 except for the Bo- 
goliubov analysis. 

The weakly interacting superfluid phase (first consid- 
ered in Ref. [28]) of the spinfull bosons is described by 
the Bogoliubov theory around a mean-field state with a 
two-component condensate wavefunction. The interme- 
diate interaction regime near the Mott transition is an- 
alyzed, as before, with a variational Gutzwiller ansatz, 
albeit with an extended local Hilbert space. However, 
unlike the spinless bosons, the ferromagnetic spin-spin 
interaction in the deep Mott phase, arising out of a 
super-exchange mechanism, is not captured by the sim- 
ple Gutzwiller ansatz. To treat this limit, we work with 
a ferromagnetic Heisenberg model with a spin-spin inter- 
action J = At 2 /U and analyze the spin-current induced 
instabilities within a spin-wave formalism. 



A. Weakly interacting Superfluid 

The weakly interacting superfluid regime admits co- 
herent mean-field spin current-carrying solutions of the 







(6) 



Such states have a spin twist of p, and carry spin cur- 
rent oc sinp • x^- between neighbors. Expanding Eq. (5) 
about this stationary state to second order in quantum 
fluctuations, b^ = biv + cj)^, gives the Hamiltonian SH = 

E k >o $ k M A where = (^jk^t-k, ^ fa-*) and 



M w = 



Un 

6- 




Un 

+ 2Un 


— e_ 



o \ 





f Un 
£ 



(7) 



Un J 



where 



e + -/i 

6 k+p/2 i e k-p/2 



with e q = — z^7 q = — 2t^ i cos^. For given U and p, 
negative eigenvalues of Mk for some k indicates LI while 
imaginary eigenvalues of aMk indicates DI, where cr = 
diag (1,-1, 1,-1). 

From here on, we will restrict ourselves to the case 
of spin currents along the diagonal of a square lattice: 
p = p (x + y) (for example p = ir represents the Neel 
state). The conditions for instabilities are summarized 
in Table I. The LI is always present for any non-zero 
pitch, while DI is always present except for the p = 
ferromagnetic state and the p = ir Neel state. 

To obtain a better understanding of the DI, we plot 
the wavevectors of the unstable modes, obtained from 
the Bogoliubov theory, as a function of the spin twist p 
in the left column of Fig. 2. Here we consider wavevec- 
tors parallel to spin current (k || p) for several values 
of U. Light gray areas indicate presence of LI but not 
DI, and dark areas indicate the presence of both LI and 
DI. The ferromagnetic state is always energetically and 
dynamically stable, as expected, while the Neel state has 
a LI but not a DI. With increasing U the region where 
the DI is present increases, i.e. more and more wavevec- 
tors become unstable. However, the region where LI is 
present is almost independent of U. 

The dispersions of the lowest collective modes (k || p) 
for three special states, the ferromagnetic state (p = 0), 
the Neel state (p = it) and the spin spiral state with a 
wavelength of 4 lattice spacings (p = 7r/4), are plotted 
in the left column of Fig. 3. For the ferromagnetic state, 
there are two low energy modes: a charge mode related 
to the U(l) symmetry breaking, which disperses linearly 
and a spin mode related to the SU(2) symmetry breaking, 
which disperses quadratically, both of which are stable 
modes. As a spin current is imposed, the charge mode 
remains stable while the spin mode develops a DI near 
k = 0, indicated by the thin red line. As we reach the 



5 



Bogoliubov Gutz wilier 




Wave Number k / n 



Figure 2: Diagrams showing regions of unstable modes for 
BEC spin-current-carrying state in optical lattice, analyzed 
with Bogoliubov theory (left column) and Gutzwiller ansatz 
(right column). The value of U/U c used is (from top to bot- 
tom) 0.1, 0.3, 0.5, 0.7. The spin twist plotted is from ferro- 
magnetic state (p = (0,0)) to Neel state (p = (7r,7r)). k is 
the wave vector of the normal mode. The shaded area in- 
dicates a LI (negative excitation energy), while dark shaded 
area indicates a DI (complex frequencies). The left column 
was previously obtained in Ref. [28]. 

Neel state, both the charge and the spin mode disperse 
linearly and are stable. Thus the DI disappears for the 
Neel state, which is stable in the weakly interacting limit. 
However, states with spin twists close to but not equal to 
7r, are unstable with the instability being seeded around 
the wave vector k = n. 



B. Gutzwiller Ansatz 

To investigate the regime near the insulator-superfluid 
transition we adopt the Gutzwiller approach of Sec. II. 
For simplicity, we restrict ourselves to the case of unity 
filling. Then there are minimally six basis states per 
site that need to be included in the local Hilbert space: 
(1°) » It) , II) , |tt) , III) , III)}, where the last three states 
have double occupancy. The local Gutzwiller wavefunc- 
tion is then parametrized in terms of ten variables (site 
indices omitted) 

\Q) = sin ^ cosy |0) +e^ 0+ ^ cos ^ cosy |t) 

+e i ^°-^ l) cos^siny ||) 



2 2 2 2 1 117 

A phase twist Ub a ) oc e* p x ) can be imposed by setting 
^0 = ^2 = P • x and other ^i/0,2 to be uniform; while 
a spin twist ((5' + ) ex e z2p x ) can be imposed by setting 
(p t = (p 4 = p . x/2 and other ^^1,4 to be uniform. Note 
that, in our parametrization, the spin-current-carrying 
states do not have any mass current, i.e. it is a state 
where the two spin species carry equal mass currents in 
the opposite direction. We find that a mass current pro- 
duces a stability diagram identical to that of the spinless 
condensate in Sec. II. From now on we will concentrate 
on the case of spin twist only. 

With a spin twist imposed on we expand the 
Hamiltonian around its stationary state and investigate 
the behavior of the fluctuation Hamiltonian. The full 
derivation is carried out in Appendix B. Here we shall 
present the results, concentrating on the case where the 
spin current is put along the diagonal: p = p^2 i Xi. For 
comparison with the Bogoliubov theory, we plot in the 
right column of Fig. 2 the wavenumber of the unstable 
wavevectors (parallel to the spin current) as a function of 
the spin twist for different interaction strengths. We find 
that, contrary to the Bogoliubov theory, the region of un- 
stable wavevectors decreases with increasing interaction 
within the Gutzwiller formalism. For example, the state 
at U/U c = 0.7 with p > tt/2 is stable in the Gutzwiller 
formalism while it shows instability within the Bogoli- 
ubov theory. The Bogoliubov theory, which is accurate in 
the weakly interacting regime, thus overestimates the dy- 
namical instability in the intermediate regime. The main 
qualitative difference, however, is in the stability of the 
Neel state (p = 7r). While the Bogoliubov theory predicts 
only a LI and no DI for this state, the Gutzwiller ansatz 
shows that the Neel state can be dynamically unstable 
in the intermediate interaction regime. 

In the middle column of Fig. 3 we plot the disper- 
sion of the low energy collective modes (with k || p) of 
the ferromagnetic state, the Neel state and a spin-spiral 
state with a period of 4 lattice spacings, for different 
interaction strengths. The ferromagnetic state has two 
gapless modes in the weakly interacting limit: a linearly 
dispersing charge mode and a quadratically dispersing 
spin mode. As interaction strength is increased towards 
the critical interaction for the superfluid-insulator tran- 
sition, /7 C , the charge mode dispersion is almost unaf- 
fected, while the spin mode dispersion flattens out. Be- 
yond the critical coupling, in the Mott phase, the charge 
mode is gapped out while the zero energy spin mode be- 
comes dispersionless. This is an artifact of our variational 
approach and we will discuss in the next section how 
this degeneracy can be lifted by considering the super- 
exchange mechanism of spin fluctuations. As soon as a 
spin current is imposed (say for the spin-spiral state), 
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Figure 3: Dispersions of excitation modes for different spin twists along diagonal p = p(x + y) where p = 0, \ and ty obtained 
with Bogoliubov theory, Gutzwiller approach, and within Heisenberg model. Solid lines are modes with real energies while the 
red dashed lines are the imaginary part of complex energies, which appears in conjugate pairs but only the positive imaginary 
part is plotted. Note the energies of collective modes from Heisenberg model are in units of 2zt 2 /U while those from Bogoliubov 
theory and Gutzwiller ansatz are in units of zt. 



the spin mode develops a DI near k = in the super- 
fluid phase. Beyond the critical coupling, the DI van- 
ishes in the Gutzwiller approach and we recover the non- 
dispersing spin mode. In the weakly interacting limit, the 
Neel state develops a dynamical instability for collective 
modes around k — n/2. This dynamical instability how- 
ever vanishes before the Mott transition point is reached. 

Comparing the results from the Gutzwiller ansatz to 
those from the Bogoliubov theory, we find that for a 
given p : the discrepancy between the two theories in- 
creases with U/t, while for a given U/t, the discrepancy 
increases with increasing p. This is understood from the 
fact that the effective Mott boundary in presence of spin 
currents is given by U = U c j p ^ and so, increasing the 
pitch of the spin-twist pushes the system closer to the 
Mott phase, where the validity of the Bogoliubov theory 
is suspect. 

Fig. 4 is the stability phase diagram of the two com- 
ponent bosons in the interaction-spin-twist plane. We 
see that any finite spin-twist leads to DI in the weakly 
interacting regime, whereas, for U > U c /2, states with 
p around tt (including the Neel state) becomes stable. 
The color scale in the plot represents the growth rate of 
the most unstable fluctuation mode in the dynamically 
unstable region. The spin 1/2 nature of the particles 
is evident in the asymmetry of the growth rate between 
p = and p = 2tt. Due to Berry's phase effects the sys- 
tem is only symmetric under 47r (and not 2tt) twist of the 
spin phase. 

As in the case of mass current in Sec. II, the valid- 
ity of our Gutzwiller approach, and hence the results 
of Fig. 4, are correct only if the fluctuation occupation 
(J2 n >o dnidni in Appendix B) is small. In the 2D case, 
we find that it is indeed small (< 0.1) for the majority 
of the stable regime, but quickly goes up near the DI 



Instability Rate 
(units of zt) 




u/u c 

Figure 4: Phase Diagram of spin current on a lattice. This is 
computed with j p = cosp which represents a spin current put 
along the diagonal of a square lattice. The shading indicates 
the growth rate of the most unstable mode in the dynamically 
unstable regime. 

boundary, which is expected as a precursor of instability. 
The most severe case happens at the DI boundary for the 
Neel state, having fluctuation occupation ~ 0.4. 



C. Heisenberg Model 

It was previously shown that the two-component Bose- 
Hubbard model (5) reduces to a ferromagnetic spin model 
in the deep Mott phase [38, 39]. Here we shall show that 
within the Gutzwiller ansatz, the ferromagnetic ordering 
is provided by the fluctuations. 

First note that since a Mott phase has = in 
our Gutzwiller ansatz, all spin twists p give the same 
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variational ground state energy. The correction to 
the ground state energy due to fluctuations is AE = 
\ S/c ( — TrM + J2 a oj a k) where M is the matrix derived 
in Appendix B and uj a k are the eigenenergies. We ex- 
pand uo a k to the first order in 1/U and find that the 
correction to the ground state energy for different spin 
twists obeys E (p) — E (0) = — (7 P — 1). This is ex- 
actly the energy difference due to different magnons in a 
ferromagnetic Heisenberg model. 

However, one should be careful in interpreting the sta- 
bility in the deep Mott regime. As noted previously [39], 
we find a non-dispersing zero mode in the Mott phase, 
which emerges because there is no energy cost to create a 
spin-flip locally. We note that this is not physical and will 
be lifted at the next order in perturbation theory. An- 
other way to understand this is that the product form we 
chose for the variational state [Eq. (8)] is unable to cap- 
ture the spin ordering in the Mott phase, because charge 
fluctuations are completely absent. To account for this, 
one can rotate the state with a suitable unitary trans- 
formation \G') = e %s |G), which amounts to a canoni- 
cal transformation H' = e~ lS He %s on the Hamiltonian 
[40, 41]. To the lowest order, H acquires the correction 
—JY^{ij) &i • Sj where J = Using this Hamiltonian 
in our analysis we find that the zero mode is indeed lifted, 
with energy 



^k(p) = Y Y ~ 7k ) ( 7 P \ ( 7k +P + 7k -p)^ (9) 

which is plotted in the right column in Fig. 3. Note 
that this dispersion is identical to the usual spin mode in 
the Heisenberg model with spin twist p. Similar to the 
Bogoliubov results, this spin mode has a LI for non-zero 
pitch and a DI for any pitch except for the ferromagnetic 
and Neel state. However a crucial difference is that the 
growth rate of the unstable modes in this case has order 

,2 

of magnitude jj, which is much smaller than that of the 
DI we find in Fig. 4. This would imply that the deep 
Mott state is at least quasi-stable, in that the instability 
time scale could be much longer than the experimental 
time scale. 



is expected to be stable [35, 44, 45] over relatively large 
time-scales, which has led to the idea that the physics of 
the Neel state may be accessed in systems which are care- 
fully prepared to be stuck in this metastable state. In the 
opposite limit of the weakly interacting superfluid phase, 
analysis using the Bogoliubov approach in Sec. Ill A also 
demonstrates that the Neel state is dynamically stable. 
This naturally leads to the question: Is the Neel state 
stable throughout the phase diagram (i.e. for all interac- 
tion strengths)? 

We use the Gutzwiller ansatz scheme to look at the 
stability of the Neel state in the intermediate interaction 
regime. The Gutzwiller approach shows that the Neel 
state is dynamically unstable for < U < 0.51. The 
U = state is technically stable, but is mostly irrele- 
vant for real experimental purposes as non-interacting 
bosons are pathological even in equilibrium (e.g. di- 
vergent compressibility) and need a finite interaction to 
form a stable superfluid. In Fig. 5, we plot the growth 
rate of the most unstable fluctuation mode (if any) of 
the Neel state, obtained via our Gutzwiller approach 
as a function of the interaction strength. The growth 
rate initially increases with the interaction strength in 
the weakly interacting limit reaching a peak at around 
U/U c ~ 0.3. It then decreases with increasing interac- 
tion and vanishes at around U/U c ~ 0.5. Thus Neel state 
physics can only be probed with dynamically generated 
states for U/U c > 0.5. We note that at U/U c = 0.55, 
the Gutzwiller ansatz predicts a condensate depletion of 
about 0.21, which shows that the approximation, which 
involves a truncated Hilbert space, captures the essential 
physics in this regime. The Gutzwiller results in the very 
weakly interacting regime, U —> 0, are, however, suspect 
as the large number fluctuations in this limit are incom- 
patible with the truncation of the Hilbert space used in 
the Gutzwiller scheme. In fact, to leading order in the 
interaction strength, the Bogoliubov theory, which pre- 
dicts a stable Neel state, is much more trustworthy than 
the Gutzwiller scheme. It would be interesting to see 
how the instability rates in the Gutzwiller approxima- 
tion change with increasing the size of the Hilbert space, 
but this much more complicated problem is beyond the 
scope of this paper. 



IV. INSTABILITIES OF THE NEEL STATE 

From the beginning of implementation of optical lat- 
tices, observing antiferromagnetically ordered states has 
been a holy grail of cold atom experiments. Although 
the original ideas involved looking for antiferromagnetic 
states with fermions, recently two component bosons 
have been proposed as an alternate medium to observe 
antiferromagnetism. In this context, there is a special in- 
terest in the observation of Neel state with a commensu- 
rate spin-ordering vector p = (tt, 7r), which is notoriously 
hard to realize as a ground state in cold atom systems 
[35, 42, 43]. In the deep Mott phase, this state is the 
highest energy state of the ferromagnetic spin model, and 



V. CONCLUSION 

In this work we have analyzed the stability of mass and 
spin current carrying two-component Bose condensates 
in optical lattices. We have approached the problem via 
Bogoliubov theory and the Gutzwiller ansatz, to handle 
respectively the weakly interacting superfluid phase and 
the regime near the Mott boundary. For small spin cur- 
rent and small interaction the two approaches agree, but 
deviate when we increase the spin current or the inter- 
action strength. In the deep Mott phase, we addressed 
the subtleties we encountered with our variational ap- 
proach, and showed that the ferromagnetic Heisenberg 
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is turned off, the system would exhibit violent disruption 
of the spin pattern if it is in a dynamically unstable state, 
as long as the inverse growth rate of the unstable modes 
are small compared to experimental timescales. Since 
the maximum growth rate of unstable modes is ~ 1.5z£, 
this growth dynamics should be observable over a wide 
range of experimental parameters. We thus hope that 
our predicted instabilities would be easily seen by future 
experiments on cold atoms. 



Figure 5: Growth rate of the most unstable fluctuation mode 
of Neel state as a function of U, the on-site interaction, ana- 
lyzed via the Gutzwiller approach. The result is expected to 
be invalid near U = because it is in superfluid phase. 



model provides an adequate description in this limit. 

For mass current carrying states, we find that the sta- 
bility phase diagram of the two component bosons qual- 
itatively follow that of the spinless bosons. The current 
carrying states are dynamically unstable beyond a critical 
value of the current, and the critical current monotoni- 
cally decreases with increase in interaction strength, van- 
ishing at the critical coupling for the superfluid-insulator 
transition. For the spin-current carrying states, we find, 
within Bogoliubov theory, that the system is unstable to 
any finite spin current in the weakly interacting limit, 
with the exception of the Neel state (with a spin twist 
of (7r,7r). The Gutzwiller approach also predicts a simi- 
lar scenario with the only difference being that the Neel 
state is also dynamically unstable in this regime. The 
Gutzwiller approach shows that the region of instability 
in the spin-current-interaction plane decreases with in- 
creasing interaction, with states around the Neel state 
(spin twist close to (ir, 7r)) being the dynamically sta- 
ble states. The Heisenberg model in the atomic limit 
also predicts dynamical instability for all current carry- 
ing states except the Neel state. Finally we stressed the 
fact that although the Neel state is stable in the two ex- 
treme limits of strong and weak interactions, it can de- 
velop instabilities for intermediate interaction strengths. 

Although energetic instabilities would be very hard to 
observe experimentally in cold atom systems, the dynam- 
ical instabilities of the current carrying states should be 
fairly easily observable as dramatic phenomena. The typ- 
ical experiment would consist of creating a spin-current 
carrying state by tuning a spatially varying artificial Zee- 
man field. Such fields with commensurate wavevectors 
have already been produced in the laboratories. By tun- 
ing the amplitude of these fields to a very large value, so 
that the Zeeman energy is the largest energy in the prob- 
lem, the spin-current carrying states can be generated as 
the ground states of the system of bosons. Once the field 
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Appendix A: Analysis for Spinless Condensate 

Here we shall give the details of calculations in Sec. II. 
The variational energy using the three-state Gutzwiller 
ansatz with phase twist p is: 



N 



(n\H\n) 



\i cos x I sin 



tzn 



-7 P sin 



1 + fT 1 sin 2 ^ 



sin x cos 



2ri) (Al) 



where 7 P = z 1 e* p-5 , which reduces to 7 P = 

cosp ^ d ~ 1 ) for a current put along an axis of a d- 
dimension cubic lattice. 

Since |0) has filling ratio v = n — sin 2 | cosx, to ensure 
commensurate filling one should find a \i such that the 
minimum of e occurs at x — f • For convenience we 
consider only the limit n ^> 1. Then with \i = 0, the 
minimum of e is attained with x = f j V = 0, an d 



1, otherwise 



(A2) 



where u = 7^7= . One can also check that this is a sta- 
tionary solution by varying e locally to leading orders 
to find Se = \ (U — 4ztnj p cos 0) sin 080. The solution 
sin# = corresponds to a Mott phase while the other 
solution corresponds to a superfluid state. 

Performing the fluctuation expansion as outlined in the 
text, one finds H = ^ Zlk ^k^k^k where, in the su- 
perfluid and Mott phase, M are respectively 
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M, 



sf 



Mr, 



cos 2 §7 + 



( 27 P — cos 2 #7 + — cos 2 #7+ — cos 6 cos |7 

— cos 2 #7+ 27 P — cos 2 #7+ — cos # cos ^7 

-cos^cos|7_ — cos^cos|7_ (27 P — 7+) cos 

y cos^cos|7_ cos^cos|7_ 

/ 2ii-7 + -7+ -7_ 7_ \ 

-7+ 2^-7+ - 7 _ 7_ 

-7- -7- 2^-7+ 7+ 

V 7- 7- 7+ 2^ - 7+ / 



2 
2 



cos # cos |7_ 
cos # cos |7_ 
cos 2 §7+ 



(2 7p -7 + )cos 2 f / 



(A3) 



(A4) 



where j± = \ (7k+p ± 7k-p)- Note that the problem re- 
duces to the one considered before [36] in the limit of 
p = 0. The spectrum is found by diagonalizing aM 
where a = diag (1, — 1, 1, — 1) which would give a spec- 
trum of the form {±|co> a k}- 



b y e, 



. cos 2 I . Then the variational energy is minimized 



>o 



and 



cos = 



r u/l P /2 ,0<ii/7p /2 <l 
I 1 , otherwise 



Appendix B: Analysis for Two- Component 
Condensate 



We take \Q) from Eq. (8) with a spin twist p imposed 
to evaluate the variational energy e = (0| H With 
some algebra it can be shown that one can first set #1,3,4 
to be 7r/2, after which 



U 



— + /i cos #2 I sin - 



(9 sin 2 (9 



"7p/2 



x (3 - cos(9 2 + 2^8^(92^ . (Bl) 
To ensure a filling ratio v = 1 one must set fi = 



where 



The six states are written in 



(3+2V2) zt' 

terms of the E: \a) — > £t |vac) satisfying the constraint 
J2 a ta^a = 1 7 where a could be the empty, 1 spin-up, 1 
spin-down, 2 spin- up, 1 spin- up + 1 spin-down, or the 2 
spin-down states. The boson creation/annihilation oper- 
ators are replaced by the pseudospin operators: 



"it 



tt t 



(B3) 
(B4) 



The unitary transformation analogous to Eq. (4) is cho- 
sen as: 



1 



"71 cos I 
1 

v/2 







ip- Xi /2 n 

-^ C0S 2 

>*P- X i/ 2 • ft 

-y2^ Sm l 








V2 





W Sm 2 

-%7T cos f 

2^2 



e ipx * 

^/2 
e* p ' x * 

V2 



1 





^V^ Sm 2 

17T C0S 2 



2^2 



"~72~ 



/4 \ 



m 
t 



V *|4i / 



(B5) 



/4 


\ 


4 




4 




4, 




4 




V4 


/ 



and we set both doi and d^ to be y 1 — X^ n >o ^lii^ni ~ 

1 — I Xl n >o dnidni because dj^ is macroscopically occu- 
pied. The validity of this expansion should be checked af- 
ter the Bogoliubov transformation to ensure consistency. 

Written in terms of d n> o and to the lowest (quadratic) 
order, H = f ]T) k (cj 4 kdj k d 4k + ^M* k ) where ^ = 



I 

d 5) -k) and M is 
an 8 x 8 matrix, whose non-zero entries are: 

Af M = M 2 , 2 = ^ (3 + 2V2) (ucosf? + 7 p/2 sin 2 0) 

- i 7+ (9 - 2^ + (3 + 2^) cos 26») 
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M 3 ,3 = M 4 ,4 = ^(3 + 2v / 2) (^ + iicos(9 + 7 p/2 sin 2 (9) = (-1 + V2 - (l + V^) cos 6>) sin ^ 

3 2 Mi 7 = M 7 1 = -M 2 8 = -M 8 2 

— -7 + cos - 

-7_ cos - (2 - >/2 + (2 + V^) cosfl) 



M 5>5 


= M 6j6 = 


M 7)7 


= M 8 , 8 = 


Mi >2 


= M 2 ,i = 


M 3)4 


= M 4 , 3 = 


Mi >3 


= M 3 ,i = 


Mi >4 


= M 4 ,i = 


Mi >5 


= M 5 ,i = 


Mi >6 


= M 6 ,i = 



1 - 2# 



—7+ (-3 + 6^ + (3 + 2^) cos26>) 



Mi, 8 


= M 8 ,i 


= -M 2)7 


= -M 7>2 






1 / 

= - 7 _ COS -(2-> 


M 3)5 


= M 5 , 3 


= -Af 4 , 6 


= -M 6i4 


M 3)6 


= M 6 , 3 


= -Af 4 , 6 


= -Af 5i4 


M 3)7 


= M 7 , 3 


= -Af 4 , 8 


= -Af 8 ,4 


M 3)8 


= M 8 , 3 


= -M 4 , 7 


= -M 7>4 


M 5>7 


= M 7i5 


= M 6>8 = 


M 8 , 6 



27+ sin - 

X - cos 2 ^ ((3 + 2^) « - 2 7+ ) M 3)5 = M 5 , 3 = -M 4l6 = -M 6l4 = ^ 

+ ^ cos 2 ^ (4 + 2^2 - (3 + 2\/2) cos6>) M 3)6 = M 6 , 3 = -Af 4 , 6 = -Af 5 , 4 = ^ sin<9 



2 p 
. cos z - 

2 



v^7+ cos 2 - M 3)8 = M 8 , 3 = -Af 4 , 7 = -M 7>4 = J= cos 2 - 

= M 2)4 = M 4>2 = -7 P /2 cos - - -7+ cos J - 

M 2)3 = M 3)2 = -7+ cos - sin 2 - 4 V 7 ,p/ 2^ 



sm( 



-M 2 , 6 = -M 6j2 

i 7 - (-1 + V2 + + cos*) sin? w J ere 7± = ^/ 1 7k+ 1 p / 1 2± 7 k 1 - p/2) ' 1 ^fS^^f^S ctM 

4 V V / / 2 where a = diag (1, —1, 1, —1, 1, —1, 1, —1) gives the spec- 

-M 2?5 = -M 5?2 trum 
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